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Abstract 

We introduce a new simplified method for computing the electron field emis- 
sion current in short carbon nanotubes using ab-initio computation in pe- 
riodic simulation cells. We computed the evolution of the wave functions 
using Time-Dependent Density Functional Theory, where we have utilized 
the Crank-Nicholson propagator. We found that in pristine carbon nan- 
otubes, the emitted charge tends to emerge mostly from electrons that are 
concentrated at the nanotube tip region. The charge beam concentrates 
into specific channel structures, showing the utility of carbon nanotubes in 
precision emission applications. 



1. Introduction 

Field emission, the escape of electrons from the surface into vacuum under 
the influence of an applied electric field, has been receiving research attention 
owing to its theoretical, as well as commercial, significance pQ. The litera- 
ture on field emission form various nanostructures is enormous. Unique field 
emission properties were discovered in carbon nanotubes owing to their high 
aspect ratio and mechanical and chemical stability, making them the best 
candidate for flat-panel displays and electron microscopy probes [21 El H]- 
Field emission from graphene was also studied [5]. Field emission was also 
investigated in nanodiamonds [6], buckyball coating EJ E] and a variety of 
nanoclusters [HUE]- Among the number of nanostructures investigated for 
field emission, carbon nanotubes proved, both theoretically and experimen- 
tally, to possess strongest field emission characteristics. 

Field emission was understood as a realization of quantum tunneling 
through the simple Fowler Nordheim theory that treats a one dimensional 
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system composed of a free electron and a trivial triangular potential barrier, 
thus without any consideration for scattering effects [12J. The escape of the 
electron through the potential barrier can be modeled using the WKB ap- 
proximation, yielding the following relationship between the extracted field 
emission current and applied electric field 

/ = 6.2 x 1(T 6 , V—- fV 2x1 ° 8x5/2 / f , 
[X + NX 1 

where \ is the thermionic work function, \i is the familiar parameter in the 
electron distribution in the Fermi-Dirac distribution function. However, at 
the nanoscale, approximating the nanotip using the sharp potential well is far 
from accurate [T3]. Thus was the need for ab-initio computational methods 
to simulate field emission at the scale of the atom. 

In recent years, field emission was simulated using a variety of ab-initio 
computational schemes, including time-dependent density functional theory 
(TDDFT) [T3l [HI EES! US] , Kubo formalism [17] , transfer matrix formalism 
[T5] and Landauer-Buttiker formalism [TJ. Here we present an adaptation of 
TDDFT. In Sec. [2] we provide a theoretical overview of our computational 
scheme, which is based on the implementation of the Octopus code [19] . We 
also indicate where we have altered the original code in order to compute the 
charge in a portion of the simulation box. 



2. Computational Details 

2.1. Ground State Calculation 

We performed Density Functional Theory computation using the Local 
Density Approximation for exchange and correlation potential by Pedrew 
and Zunger [20J. States are expanded in terms of plane waves, 

^(r) = ^^(G)e iGT 
G 



^n(G) = i J dr^ n (r)e 4G ' r , 



where V is the super-cell volume, and ip n (r) and ip n (G) are related by a 
three-dimensional Fourier transform. Atomic orbitals are substituted with 
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Troullier-Martins pseudopotentials. The density matrix is mixed according 
to the Broyden mixing scheme [21J. We perform the ground state calculation 
in the absence of the external electric field, and turn on the field during the 
time-dependent calculation. 

2.2. Time-Dependent Density Functional Theory 

The solution of the time-dependent Kohn-Sham equations 

V 2 



has the following exact expression: 



+ v K s(r,t) 



^ . (T) = Texp | -ijf drH(r) j ip 



,(°) 



where Texp is the time-ordered exponential, which is a short-hand for: 

ut) = |e jf dTl ' ' ' jf dTnH (Tl) ' ' ' ^ (Tn) } ^ (0) 



. n=0 



If the Hamiltonian commutes with itself at different times, we can drop 
the time-ordering product, and leave a simple exponential. If the Hamilto- 
nian is time- independent, which makes it trivially self commuting, the solu- 
tion is simply written as: 

ipj(T) = exp {—iTH} ipfK 

Unfortunately, this is not the case for TDDFT when the system is exposed 
to external time-dependent perturbations like electric and magnetic fields or 
pulsed lasers. But even without an external time-dependency, there remains 
the intrinsic time-dependency of the Kohn-Sham Hamiltonian, which is built 
" self-consistently" from the varying electronic density. 

The first step to tackle this problem is to split the propagation of the long 
interval [0, T] into iV smaller steps by utilizing the group-theoretic property 

U(T,t) = U(T,t')U(t',t) 
of the time evolution operator. This yields the following time discretization: 
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N-l 



U(T,0)= Y[u(U + At,U), 



i=0 



where t = 0, t N = T, At = T/N. So at each time step we are dealing with 
the problem of performing the short-time propagation: 

ipj(t + At) = U(t + At,t)ipj(t) = Texp l-i J drH(r)jiJ) S (t). 

2.3. Propagation Method 

There is a wide range of methods for computing the propagator for the 
time-dependent Kohn-Sham equations, and research in this area continues 
to grasp more attention [22]. These propagators solve the problem of ap- 
proximating the orbitals ipj(t + At) by the knowledge of ipj{r) and H{r) for 
< t < t. Some methods require the knowledge of the Hamiltonian at some 
points r in time between t and t + At. We choose to utilize the Crank- 
Nicholson method, also known as the implicit midpoint rule. This method 
is the most famous propagation method used to compute the time-evolution 
of the Schrodinger equation, and is based on the following implicit midpoint 
rule [231 I2H 122], which is the average of the forward Euler and backward 
Euler integration methods: 

l ^7 = H(t n+1/2 )- {lp n+1 + Ipn) , 

where 

tn+i/2 = 1/2 + t n ) ,t n = nAt. 

Then, by straightforward algebra we obtain the following: 

, 1 - iAt\H(t n+1/2 ) 

1 + iAt^H{t n+l / 2 ) 

which is also known as the Cayley approximation of the exponential to the 
second order of At. To the third order of At, the Cayley approximation is 
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1 - iAt\H{t n+1/2 ) 

1 + iAt\H{t n+1/2 ) + U[m [ 
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The problem of evaluating the evolution operator becomes that of solving 
the following linear matrix equation: 



At , 
I + i—H{t n+l/2 ) 



At . . 
I - i^H(t n+1/2 ) 



This method has two important properties: it preserves the norm of the 
wave function at all times (unitarity property), and exchanging n <H- n + 1 
gives the same numerical results (time reversibility property). However, there 
is a tendency to prefer the other method, which is the Split-Operator method 
where the exponential of the Hamiltonian is split into exponentials of its 
composite terms (the kinetic and potential energies) [22] • Nevertheless, good 
accuracy can then be achieved by the Crank-Nicholson method only for very 
small time steps (such as the time step used here, which is 0.00242 fs, as the 
computational error is proportional to the cube of the time step at each step 
of the algorithm [251 [23]). A shortcoming in the Octopus code we are using 
is that the code uses the retarded hamiltonian H(t n ) in place of H{t n+ i/ 2 ) in 



At , 
I + i—H{t n+l/2 ) 



ip{n + l) 



At . . 



which is essentially the predictor step's equation: 



I + i^H(t r 



•At „. 
i — H(t r 



The consequence of using the retarded Hamiltonian is the continuous decrease 
in energy every time-step because we will be using a retarded potential that 
arises from the retarded Hamiltonian. 



2.4- Current Calculation 

The amount of the electron remaining in the emitter region at time t for 
a particular wave function is given by 



Qn(t)= ^ J J WdXL, 

and in terms of the life time r n of state n, 

Q n (t) = e~^, 



(1) 



(2) 
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which implies linear behavior of Q n (t) in the short time interval. From this 
formula, the current generated by a wave function is given by 



dQnit) 



1 



(3) 
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= e 



dt 



r, 



n 



in the short time range. The total current is 



n 



dQnif) 

dt 



In order to compute the charge remaining in the lower portion of the sim- 
ulation box, we integrate that lower portion each time-step for each wave 
function. At the beginning of the simulation, the charge is typically 1 (which 
works as a check for normalization consistency of the propagator used). 
Charge starts to decrease as time elapses, until it reaches a minimum point 
after which it starts increasing (due to reflection of the wave back from the 
upper surface of the simulation box). The integration is performed by sim- 
ply summing the product of unit volume boxes of the mesh used in Octpus 
within the region concerned. 

3. Application to short pristine carbon nanotubes 

Field emission from carbon nanotubes has been extensively studied from 
theoretical and experimental perspectives, owing to their unique emission 
properties. Theoretical ab-initio computation of field emission from carbon 
nanotubes has investigated a number of considerations that we would like 
to highlight from the literature. Han et al. [T5] stated that using a Heav- 
iside external field results in unwanted oscillations in the charge evolution, 
and thus they prefered to use a constant field (converging the ground state 
with an external electric field) and then correct the wave functions by us- 
ing a prescribed method. Screening, the damping of electric fields caused 
by the presence of mobile charge carriers (neighboring carbon nanotubes in 
a nanotube bundle), were studied by Chen et al. [25]. They discovered an 
interplay between the spacing between CNTs and their electronic structure. 
When the array spacing between neighboring nanotubes is three times the 
nanotube length, the applied external field is strongly screened. Most simu- 
lations in the literature were performed on 2D periodic systems (an infinite 
slab structure) with infinitely many neighboring nanotubes. Charging the 
nanotube (adding extra electrons) was investigated by Lou et al. [TT] to 
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Figure 1: Charge evolution for the highest energy level upon applying a field of 0.8V/A. 
Charge is computed as the integration of the modulus of the respective wave function 
within the lower region of the nanotube, and units are arbitrary. 

simulate the realistic situation, who have also performed the simulation on 
isolated pristine nanotubes with 50-60 carbon atoms. Although the results 
reported from each ab-initio simulation methods would give different results 
for the extracted current, they agree that field emission from carbon nan- 
otubes arise mainly from states right below the Fermi level, not from the 
metallic continuum as in usual metallic emitters [27J. 

To demonstrate our method, we studied field emission from a short (5,5) 
carbon nanotube, which is a metallic conductor. Our nanotube is composed 
of 70 carbon atoms, where the bottom carbon atoms are passivated with 
hydrogen atoms in order to avoid the presence of unnecessary dangling bonds. 
The hydrogen-passivated nanotube is 8.6A long, with a diameter of 6.75A. 
The simulation box size is 30A x 30A x 30A. The tip of the nanotube lies at 
the center of the simulation box, and we calculate the charge in the lower 63% 
of the simulation box (which includes the lower 50% plus an extra volume 
past the nanotube tip). 

We adopted the time step of 0.1 au = 0.00242 fs, the same one used in 
|13j . We also experimented with a smaller time step (0.000242 fs) and we 
obtained the same results, which shows that this time step is optimal for the 
Crank-Nicholson algorithm. In Tab. [T] below, we show comparison between 
our results and those reported in Ref. [13j and Ref. [15] (the cited quantities 
are extrapolated from graphs included in the references). 



Our simulation results in a periodic simulation box are close to the results 



Table 1: Extracted current (//A) from pristine carbon nanotube against an applied electric 
field (V/A). Our results are almost close to those reported in [IS] when the applied field 
is in the range 0.2-0.8V/A. For larger values, reported results for the current differ from 
the result reported in |13) due to the difference in length of the nanotube used. 
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reported by Han et al. (13] for the applied field in the range 0.2V/A to 
0.6V/A, and are very close to the result reported by Mayer et al. [28J at a 
field of 0.25V/A (note that Han et al. have reported different results in two 
different publications, |T3j and [15] as is indicated in the Tab. [I] above). At 
an applied field of 1.0V/A, both results reported by Han et al. are slightly 
higher than ours, which is due to the difference in size between the nanotubes 
used; that is, shorter nanotubes generally yield lower current than longer 
nanotubes [T4] . 

In Fig. [2j below, we present the evolution of the cross section of the 
highest energy wave functions in a pristine nanotube under a external field 
of 0.8V/A. 

We notice in Fig. [2] the evolution of charge of the wave function as the 
base of shape (A) extends forward. We also note that by 15 branch 
emerges (B) from the shape, which vanished completely by 30 the 
wave function becomes more concentrated into the shape (A) . Also note that 
the relative height of shape (A) is less at time 30 a.u. than that at time 5 a.u., 
which is due to the fact that a portion of the charge of that wave function 
flows into the vacuum region. We can gain a better insight by noticing the 
evolution of the terminal curves at the end of the simulation box below. 
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Figure 2: Evolution of the wave function into the vacuum region: highest energy wave 
function. Units are arbitrary. 
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Figure 3: Evolution of the wave function into the vacuum region: highest energy wave 
function, close to the upper surface of the simulation box. The horizontal axis is the 
intersection between the plane crossing the nanotube vertically and the upper face of the 
simulation box. The vertical axis shows the charge computed in Eq. [2] above, in arbitrary 
units. 

Fig. [3] shows clearly that the wave functions are concentrated into cer- 
tain "channels" as time elapses, which indicates that the emitted current 
is concentrated towards the anode, instead of being dispersed into space. 
Studying such behaviors can be instrumental when the exact shape of the 
electron beam is concerned (when producing very high precision beams in 
atomic force microscopes is required). It is also of great importance whether 
we can influence the shape of the beam by adding certain dopants, or by 
performing structural modifications. This issue will be treated in a later 
publication. 

In such a short nanotube, dangling bonds at the terminal carbon atoms 
(opposite to the tip) play a significant role. We observed that field emis- 
sion in the non-hydrogen-passivated isolated carbon nanotube is almost 50% 
higher than that of the hydrogen-passivated one. This is in line with the 
idea established in [TBj, that dangling bonds (arising from vacancy defects 
and non-passivation with hydrogen at the terminal bonds) are the prime 
contributors to the field emission current. 

4. Conclusions 

Using time-dependent density functional theory, we studied the mecha- 
nism of field emission from a short carbon nanotube. We utilized the classi- 
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cal Crank-Nicholson method to compute the wave function propagator, and 
used a small time step in order to reduce the computational error. Dangling 
bonds tend to considerably increase the field emission current in nanotubes; 
7r bonds do not contribute to field emission as much as a bonds. We found 
that the higher electronic density at the nanotube tip (which takes place due 
to the pentagonal structure, and which occupy higher molecular orbitals) 
are the prime contributors to field emission, unlike the case of metallic tips 
where emission is produced by the continuum of electrons in the metal [27]. 
We reproduced the results reported by Han et al. [13] by using periodic 
boundary conditions. The charge cloud emerging from the tip appears to 
concentrate into narrower channels as the cloud approaches the opposite end 
of the simulation box. 
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